arXiv:1506.02992v2 [cond-mat.soft] 15 Jun 2015 


Macroscopic model with anisotropy based on micro-macro 

informations 

N. Kumar, S. Luding, and V. Magnanimo 

Multi Scale Mechanics (MSM), Faculty of Engineering Technology, MESA+, 
University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands 

June 16, 2015 

Abstract 

Physical experiments can characterize the elastic response of granular materials in terms 
of macroscopic state-variables, namely volume (packing) fraction and stress, while the mi¬ 
crostructure is not accessible and thus neglected. Here, by means of numerical simulations, 
we analyze dense, frictionless, granular assemblies with the final goal to relate the elastic 
moduli to the fabric state, i.e., to micro-structural averaged contact network features as 
contact number density and anisotropy. 

The particle samples are first isotropically compressed and later quasi-statically sheared 
under constant volume (undrained conditions). From various static, relaxed configurations 
at different shear strains, now infinitesimal strain steps are applied to “measure” the effective 
elastic response; we quantify the strain needed so that plasticity in the sample develops as 
soon as contact and structure rearrangements happen. Because of the anisotropy induced by 
shear, volumetric and deviatoric stresses and strains are cross-coupled via a single anisotropy 
modulus, which is proportional to the product of deviatoric fabric and bulk modulus (i.e. 
the isotropic fabric). Interestingly, the shear modulus of the material depends also on the 
actual stress state, along with the contact configuration anisotropy. 

Finally, a constitutive model based on incremental evolution equations for stress and 
fabric is introduced. By using the previously measured dependence of the stiffness tensor 
(elastic moduli) on the microstructure, the theory is able to predict with good agreement 
the evolution of pressure, shear stress and deviatoric fabric (anisotropy) for an independent 
undrained cyclic shear test, including the response to reversal of strain. 

Keywords: Anisotropy, shear, calibration, prediction, PARDEM. 


1 Introduction 

Granular materials behave differently from usual solids or fluids and show peculiar mechanical 
properties like dilatancy, history dependence, ratcheting and anisotropy inni [Till EH cni inni sni 
[THlleollZQlEHllHZ!. The behavior of these materials is highly non-linear and involves plasticity 
also for very small strain due to rearrangements of the elementary particles |22l [T^ Uj. The 
concept of an initial purely elastic regime (small strain) for granular assemblies is an issue still 
under debate in the mechanical and geotechnical communities. On the other hand, approaches 
that neglect the effect of elastic stored energy, i.e., all the work done by the internal forces 
is dissipated, are also questionable. Features visible in experiments, like wave propagation, 
can hardly be described without considering an elastic regime. In a general picture, both the 
deformations at contact and the irrecoverable rearrangements of the grains sum up to the total 
strain. The former represents the elastic, reversible contribution to the behavior of the material. 
That is, for very small strain the response of a finite granular system in static equilibrium can 
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be assumed to be linearly elastic [591 Eol EHl 1121, as long as no irreversible rearrangements take 
place. 

Despite these arguments and the long-standing debate, basic features of the physics of gran¬ 
ular elasticity are currently unresolved, like the determination of a proper set of state variables 
to describe the effective moduli. Physical experiments carried out on sand and glass beads show 
that wave propagation in the aggregate depends upon the stress state and the volume fraction 
[laiMiisniEaEaiis]. Recent works [Ml 113 [Ml EZl [I] show that along with the macroscopic 
properties (stress and volume fraction) [Ml EHl ES] , also the fabric [6ll 071 [73 EZl [TO] plays a 
crucial role, as it characterizes, on average, the geometric arrangement of contacts. Due to prepa¬ 
ration and loading path, the microstructure of granular aggregates is often far from isotropic 
and this is at the origin of interesting features in those materials. The mechanical behavior 
of anisotropic soils is a topic of current interest for both experimental and theoretical investi¬ 
gations. As one example, extensive experimental work of anisotropy has been carried out on 
laboratory-prepared (by careful ‘raining’ or bedding) sand specimens [STIEI]- These and other 
studies show that the sample deformation characteristics depends highly on the orientation of 
the bedding plane with respect to the principal stress and strain axes [13 ED EH [83 Ea @3] On 
the other hand, when the material is sheared, anisotropy in the contact network develops, as 
related to the opening and closing of contacts, restructuring, and the creation and destruction 
of force-chains, affecting the material response. [3 [771 EH @3 E3. 

Most standard constitutive models, involving elasticity and/or plasticity have been applied 
to describe the incremental behavior of (an)isotropic granular solids - sometimes with success, 
but typically only in a limited range of parameters. In the majority of the models, the stress 
increment is related to the actual stress state of the granular system and its density. This 
is the case for hypoplasticity |23l I36| . where a single non-linear tensorial equation relates the 
Jaumann stress-rate with strain-rate and stress tensors. Only few theories after the pioneer work 
by Cowin [TO], consider explicitly the influence of the micro-mechanic structure on the elastic 
stiffness, plastic flow-rule or noncoaxiality of stress and strain, see [73 [73 E3 03 EH B E] and 
references therein. The evolution of microstructure due to deformations is an essential part of 
a constitutive model for granular matter because it stores the information how different paths 
have affected the mechanical state of the system. In this sense, fabric is a tensorial history 
variable. When included in the formulation, the effect of structure is often described by a 
fixed fabric tensor normal to the bedding plane of deposited sands [73 03 EH @3. Recently 
Li & Dafalias 03 have proposed a new framework (rather than a specific constitutive model) 
by reconsidering the classical steady state theory by Roscoe et al. [63, with a fabric tensor 
evolving towards a properly defined steady state value. This is supported by experimental [M] 
and extensive numerical works [281 (301 071 IFTl [77] • In a similar fashion, the anisotropy model 
proposed in 0301] postulates the split of isotropic and deviatoric stress, strain and fabric and 
includes the microstructure as a variable, whose behavior is described by an evolution equation 
independent of stress. Refs. [Ml [33 predicts uniaxial simulation results under this assumption 
(independent evolution of stress and structure), where the simplified model well captures the 
qualitative behavior. 

In this work we use the Discrete Element Method (DEM) to study granular assemblies made 
of polydisperse frictionless particles and focus on their elastic behavior. By isolating elasticity 
we aim to distinguish the kinematics at the microscale that lead to either macroscopic elasticity 
or plasticity. We analyze the role of microstructure, stress state and volume fraction on the 
evolution of the elastic moduli, with the goal to characterize them in terms of a unique, limited 
set of variables. In order to calculate the stiffness tensor, we apply small-strain probes to various 
equilibrium states along a volume conserving (undrained) shear deformation path. In the case 
of a finite assembly of particles, in simulations, an elastic regime can always be detected and 
the elastic stiffnesses can be measured by means of an actual, very small, strain perturbation 
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[50] , The purpose is to improve the understanding of elasticity in particle systems and to guide 
further developments for new constitutive models. As an example, the relation between moduli 
and fabric here is used in the anisotropic constitutive model, as proposed in Ha El], to predict 
the macroscopic behavior during a more general deformation path, involving also strain reversal. 

This paper is organized as follows: The simulation method and parameters used and the 
averaging definitions for scalar and tensorial quantities are given in section [21 The preparation 
test procedures, and the results from the deviatoric simulation are explained in section El Section 
m is devoted to the measurement of elastic moduli by means of small isotropic and deviatoric 
perturbations. There we present the evolution of the moduli with strain and link them to 
fabric and stress. Finally, section O is devoted to theory, where we relate the evolution of the 
microstructural anisotropy to that of stress and strain, as proposed in Refs. HSl El]. This 
displays the predictive quality of the model, calibrated only for isochoric, uni-directional shear, 
when applied to an independent, cyclic shear test. 

2 Numerical simulation 

The Discrete Element Method (DEM) [21 ETl EQ] helps to study the deformation behavior of 
particle systems. At the basis of DEM are laws that relate the interaction force to the overlap 
(relative deformation) of two particles. Neglecting tangential forces, if all normal forces fj acting 
on particle i, from all sources, are known, the problem is reduced to the integration of Newton’s 
equations of motion for the translational degrees of freedom: 

^ (rriiVi) = fi -k rriig, (1) 

with the mass m* of particle i, its position r^, velocity Vj (= r-j) and the resultant force fj = fj'’ 
acting on it due to contacts with other particles or with the walls, and the acceleration due to 
gravity, g (which is neglected in this study). The force on particle i, from particle j, at contact 
c, has normal and tangential components, but the latter are disregarded in this study to focus 
on frictionless packings. 

Eor the sake of simplicity, the linear visco-elastic contact model for the normal component 
of force is used, 

r = k6 + 7 < 5 , ( 2 ) 

where k is the spring stiffness, 7 is the contact viscosity parameter, 5 = {di + dj) /2 — (r* — r^) • 
h is the overlap between two interacting species i and j with diameters di and dj, h = 
(rj — Fj) / |(rj — rj)| and 5 is the relative velocity in the normal direction. In order to reduce 
dynamical effects and shorten relaxation times, an artificial viscous background dissipation force 
ffe = —7feVi proportional to the moving velocity Vj of particle i is added, resembling the damping 
due to a background medium, as e.g. a fluid. 

The standard simulation parameters are, N = 9261(= 21^) particles with average radius 
(r) = 1 [mm], density p = 2000 [kg/m^j, elastic stiffness k = 10® [kg/s^j, particle damping 
coefficient 7 = 1 [kg/s], background dissipation 7^ = 0.1 [kg/s]. Note that the polydispersity of 
the system is quantified by the width {w = rjnax/''’min = 3) of a uniform size distribution [M] , 
where rmax and are the radii of the biggest and smallest particles respectively. 

The average time time scale is determined when two averaged size particle (with ravg = 
(r) = 1) with mass mavg = p (dvrravg/S) = 8.377 [^g] interact, and is given as tc,avg = 

^/-y/^/^-avg “ (7/ (2"iavg))^ =0.6431 [ps], where is the reduced mass, with 

restitution coefficient 

Cavg = exp (2mavg)) = 0.926. The fastest response time scale in the system is deter¬ 
mined when two smallest particle with mass msmaii = P ~ 1.047 [^g] interact, and is 
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given as tc,small = Tr/y/fc/m^^^n - ( 7 / = 0.2279 [/xs], where = msmaii/2 is the 

reduced mass, with restitution coefficient esmaii = exp {—'jtc ,small/ — 0.804. 

2.1 Coordination number and fraction of rattlers 

In order to link the macroscopic load carried by the sample with the active microscopic contact 
network, all particles that do not contribute to the force network are excluded. Frictionless 
particles with less than 4 contacts are thus ‘rattlers’, since they cannot be mechanically stable 
and hence do not contribute to the contact or force network [30112^ |39] . The classical definition 
of coordination number is C = M/N, where M is the total number of contacts and N = 9261 is 
the total number of particles. The corrected coordination number is C* = where, M 4 is 

the total number of contacts of the A ^4 particles with at least 4 contacts. Moreover, we introduce 
here the reduced number of contacts where contacts related to rattlers are excluded twice, 
as they do not contribute to the stability of both the rattler and the particle in contact with 
it. Hence, = M 4 — Mi — M 2 — M 3 = M — 2 (Mi + M 2 + M 3 ), where Mi, M 2 and M 3 are 
total number of contacts of particles with only 1, 2 and 3 contacts respectively. This leads to 
a modification in the definition of the corrected coordination number is C* = M|’/A^ 4 . The 
fraction of rattlers is (j)r = (N — N 4 ) /N, hence, C = C* (1 — (j)r)- The total volume of particles 
is ^ ~ 4:TrN{r^), where (r^)/3 is the third moment of the size distribution [241139j and 

volume fraction is defined as = (1/H) Ylv=i where V is the volume of the periodic system. 

2.2 Macroscopic (tensorial) quantities 

Here, we focus on defining averaged tensorial macroscopic quantities - including strain-, stress- 
and fabric (structure) tensors - that provide information about the state of the packing and 
reveal interesting bulk features. 

By speaking about the strain-rate tensor E, we refer to the external strain that we apply to 
the sample. The isotropic part of the infinitesimal strain tensor Sy [301 IMl 139] is defined as: 

Ssy = -e^dt = _^ _ltr(hE) = -■^tr(E)dt, (3) 

where eaa= Eoadt with aa = xx, yy and 22 ; are the diagonal components of the tensor in the 
Cartesian x — y — z reference system. The trace integral of 3ev is denoted as the volumetric 
strain e^, the true or logarithmic strain, i.e., the volume change of the system, relative to the 
initial reference volume, Vq. 

On the other hand, from DEM simulations, one can measure the ‘static’ stress in the system 
[ 2 ] as 

a = (l/H)^F®r, (4) 

cev 

average over all the contacts in the volume V of the dyadic products between the contact force 
f'’ and the branch vector F, where the contribution of the kinetic fluctuation energy has been 
neglected HZlEo]. The isotropic component of the stress is the pressure P = tr((T)/3. 

In order to characterize the geometry/structure of the static aggregate at microscopic level, 
we will measure the fabric tensor, defined as 

^ vev ceT 

where V'^ is the volume relative to particle V, which lies inside the averaging volume H, and n'^ 
is the normal unit branch-vector pointing from center of particle V to contact c |47l [40l [ 86 ] . We 
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want to highlight that a different, convention for the fabric tensor involves only the orientation 
of contacts as follows [MlEIllST!: 

F° = — V (2) n" ( 6 ) 

Nc ^ 

^ ceN^ 


where Nc is the total number of contacts in the system. An approximated relationship between 
Eqs. dS]) and Q can be derived as: 



(7) 


with tr(F‘’) = 1. This relation is exactly equal for monodisperse assemblies but largely deviates 
for assemblies with high polydispersity (see further discussion in section [3]). The difference also 
becomes more significant when the jamming volume fraction 1S3 US] is approached. In the 
following, when not explicitly stated, we will refer to Eq. ([5]), since we combine the effects of 
volume fraction and number/orientation of contacts, both relevant quantities when the elastic 
moduli are considered |24j . 

In a large volume with a given distribution of particle radii, the relation between the isotropic 
fabric, i.e., the trace of F, is proportional to the volume fraction v and the coordination number 

C Refs. [Ml 121 EH] as 

= tr(F) = g^uC = gsuC* {1 - cj)r), ( 8 ) 


where C, C* and (pr have been introduced in previous section [ 2 T] and g^ ~ 1.22 for polydispersity 
w = 3, being only a weighted, non-dimensional moments of the size distribution |2a[ZIlEH|. 


2.3 Isotropic and Deviatoric parts 

We choose here to describe each symmetric second order tensor Q, in terms of its isotropic part 
(first invariant) and the second 


J2 = \ [{Q?f + + iQ^f] 

and third 

J3 = det(Q'') = Qf Qi'Qf, 

invariants of the deviator, with QiiQ^ and eigenvalues of the deviatoric tensor = 
Q — (tr(Q)/3)I. We use the following definition (of the Euclidean or Frobenious norm) to 
quantify with a single scalar the magnitude of the deviatoric part (JOllM] of Q: 


Qdev = Fsgn (Q) v^2^ = Fsgn (Q) 


{Qxx — Qyy)^ + {Qyy ~ Qzzf' + {Qzz — Qzz)^ + 6 {Q 


xy 


+ Q 2 


yz 


+ Qlx) 


(9) 

where Qxx, Qyy and Qzz are its diagonal, and Qxy, Qyz and Qzx its off-diagonal components 
and the deviators edev; CTdev and Fdev refer to strain E, stress a and fabric F, respectively. 
Fsgn(Q) is the sign function that relates the tensorial quantity to be measured, Q, with the 
reference-tensor that describes the (strain- or stress-controlled) path applied to the sample, Hq: 


Fsgn(Q) = sgn(Ho : Q). 


For a given, complex deformation path, the reference tensor Hq must be chosen in a convenient 
way, in order to take into account both the actual loading path and/or the previous deformation 
history of the sample. In the special case of undrained shear test, as introduced later in section 
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[3l we use as reference Hq = —E = (—1 ,1,0), where only the diagonal values are given, so that 
Fsgn simplifies to 

Fsgn (Q) = sgn {Qyy - , 

with X—wall expanding, y—wall compressing and z—wall non-mobile [39]. We want to point out 
here that, during a deformation, the response of stress <t and fabric F is opposite in sign to 
applied strain rate E. Unless mentioned explicitly, we will be using a sign convention for strain 
(isotropic = —(l/3)tr((iE) and deviatoric fedev = “<^-Edev)) such that consistently a positive 
strain increment leads to a positive stress and fabric response. 

Finally we note that in this work, we will use k* = k/ (2(r)) to non-dimensionalize pressure 
P and deviatoric stress (Jdev to give P* and cr^g^, respectively, and will be referring to deviatoric 
stress as shear stress. 0 

3 Volume conserving (undrained) biaxial shear test 

In this section, we first describe the sample preparation procedure and then the details of the 
numerical shear test. 

The initial configuration is such that spherical particles are randomly generated, with low 
volume fraction and rather large random velocities in a periodic 3D box, such that they have 
sufficient space and time to exchange places and to randomize themselves. This granular gas is 
then compressed isotropically, to a target volume fraction uq = 0.640, sightly below the isotropic 
jamming volume fraction [521 [3(11 (39] Vc ~ 0.658 and then relaxed to allow the particles to 
fully dissipate their potential energy [301139] . H The relaxed state is then compressed (loading) 
isotropically from uq to a higher volume fraction ofu = 0.82, and de-compressed back (unloading) 
to uq [30l [39] . 

The preparation procedure, as described above provides many different initial configurations 
with volume fractions I'i, each one in mechanical equilibrium. Starting from various r'j chosen 
from the unloading branch [301 I39| . the samples are then sheared keeping the total volume 
constant, that is with a strain-rate tensor 


E — £dev 


-10 0 

0 0 

0 0 0 


( 10 ) 


where edev = 28.39 [/rs“^] is the strain-rate (compression > 0) amplitude applied to the moving 
X— and y—walls, while the third wall is stationary. Our shear test, where the total volume 
is conserved during deformation, resembles the undrained test typical in geotechnical practice 
m- The chosen deviatoric path is on the one hand similar to the pure shear situation, and on 
the other hand allows for simulation of the biaxial element test [561 [55] (with two walls static, 
while four walls are moving, in contrast to the more difficult isotropic compression, where all 
the six walls are moving). Pure shear is here used to identify constant volume deviatoric loading 

^It is important to point out that the rattlers are excluded in defining the (corrected) coordination number 
C*. However dynamic rattler particles with 1 < Mp < 3 contacts are included in the definitions of fabric and 
stress. We verified that during shear deformation, the maximum contribution in deviatoric stress due to rattlers 
is 0.03%, while in the case of deviatoric fabric the contribution rises to 0.5%. This is not surprising since only 
contacting particles contribute to the definitions of both stress and fabric and dynamic rattlers have a smaller 
weight for stress than for fabric, see Eq. (0. Note also that the number of rattlers decreases with increasing size 
of the particles m- 

^Note that the jamming volume fraction is given for a uniform radius distribution for polydispersity w = 3. 
The results will be different if the distribution is different, e.g., when uniform surface or volume distributions are 
used. See Ref. m for a detailed discussion on the evolution of jamming volume fractions with polydispersity for 
a uniform radius distribution. 
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with principal strain axis keeping the same orientation as the geometry (cuboidal) of the system 
for the whole experiment. In this case, there is no rotation (vorticity) of the principal strain 
(rate) axis and no distortion/rotation of the sample due to shear deformation. Different types of 
volume conserving deviatoric deformations can be applied to shear the system, but very similar 
behavior has been observed [30], in terms of shear stress. 


3.1 Evolution of stress 


The evolution of non-dimensional pressure P* with deviatoric strain e^ev is presented in Fig. 


1(a) during undrained shear tests for some exemplary volume fraction. For frictionless systems 
analyzed here, only a slight variation of the pressure is observed at the beginning of the test, 
due to the development of anisotropy in the sample, after which P* remains constant]! Both 
the (small) initial pressure change and the final saturation value vary with the vicinity of i' to 
the jamming volume fraction Vc- Interestingly, depending on the volume fraction, some of the 
samples show increase of the pressure (dilatancy) with respect to the initial value and some other 
decrease (compactancy), as shown in Fig. [TJ This supports the idea of a certain threshold value 


= 0.79, as shown in Fig. 2(a) where the pressure of the system changes behavior, similarly 
to the switch between volumetric dilation and contraction visible in triaxial tests. 

The evolution of the (non-dimensional) shear stress during shear, as function of the 
deviatoric strain Sdev) is shown in Fig. 1(b), for the same simulations as in Fig. 1(a) The stress 


grows with applied strain until an asymptote (of maximum stress anisotropy) is reached where 
it remains fairly constant - with slight fluctuations around the maximum [T0|. The growth 


rate and the asymptote of (r^g.^, both increase with u. 



Figure 1: Evolution of non-dimensional (a) pressure P* and (b) shear stress cr^g.^^ along the main 
strain path for the pure shear deformation mode for five different volume fractions, as given in 
the inset. 


3.2 Evolution of fabric 

Complementary to stress, in this subsection we study the evolution of the microstructure in 


the sample during the volume conserving shear test. Fig. 3(a) shows that the isotropic fabric 
Fv behaves in a very similar fashion as P*, with a slight increase/decrease at the beginning, 


®We observe a much more pronounced change in pressure when friction is included in the calculation, in 
agreement with other studies, see e.g. |28| . These data are not shown here and are subject of ongoing research. 
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Figure 2: Difference between the final and initial values in (a) non-dimensional pressure P* and 
(b) isotropic fabric for the pure shear deformation mode for different volume fractions. Red 
represents the change in bulk modulus, as derived in section 14.31 Dashed lines in the plots 
represent the crossover when these quantities change sign. 


followed by saturation stage, whose value increases continuously with v. Fig. |2(b)| shows that 
the difference between the initial value of and its saturation value, changes sign when a 
certain volume fraction, = 0.755, is reached. Note that ^ that further confirms the 
independent evolution of F and cr. 

From Eq. ([8]) F^ is proportional to the product of volume fraction ly, that remains unchanged 
during deviatoric deformations, and coordination number C, that varies only slightly for sheared 
frictionless systems [30]. Note that as C = C* {1 — 4>r), knowing the (empirical) relations of C* 
and (pr with volume fraction, as presented in Refs. [301139| . we can fully describe the isotropic 
fabric state. In this study, we assume F^ to stay constant during the shear test. This assumption 
will be used later in section[5]for the prediction of a cyclic shear test. However, the small changes 
in Ev or P* can be associated to a (small) change in the jamming volume fraction |41j . 

The evolution of the deviatoric fabric, E^ev! as function of the deviatoric strain is shown in 
Fig. |3(b)| during shear for five different volume fractions. It builds up from different random 
small initial values (due to the initial anisotropy in the sample that develops during preparation) 
and reaches different maxima. The deviatoric fabric builds up faster at lower volume fractions 
but the maximal values are higher for smaller volume fractions, qualitatively opposite to the 
evolution of [TU] . As mentioned in section 12.21 the validity of Eq. (ED, that relates the two 
different definitions of fabric depends on polydispersity. In order to check the relation, in Eig. 
[4] the evolution of the three eigenvalues of the fabric tensor is plotted, for both definitions, Eqs. 
dSD and Q, during the volume conserving shear test, for three different values of polydispersity 
w =1, 2 and 3. Eor all polydispersities, the chosen volume fraction u = 0.685 is close to the 
jamming points, that slightly varies with w [39|. The difference between the definitions of 
fabric becomes higher for higher polydispersity tc = 3, as in Eq. (|6D the contribution of each 
particle is weighted to its surface area, whereas in Eq. ([5D it is weighted by the volume. Only 
for the monodisperse case, the relation is exact, as can be seen in Eig. 4(a), The differences 
are considerable for w = 2 and w = 3, for both compressive and tensile direction, while the 
non-mobile direction is not affected. Note that the difference of the two fabrics will be smaller 
for denser systems. 




















Figure 3: Evolution of (a) isotropic fabric and (b) deviatoric fabric F^ev along the main 
strain path for the pure shear deformation mode for five different volume fractions, as given in 
the inset. 
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Figure 4: Evolution of the eigen-values of the fabric tensors (directions shown in the inset), 
during shear deformation at volume fraction u = 0.685, for the fabric dehnition defined in Eq. 
(j6|) (smaller symbols) and the relation presented in Eq. ([7]) (large symbols), for three cases of 
polydispersity (a) rc = 1, i.e., monodisperse (h) w = 2 and {c) w = 3 (present work). 


4 Elastic moduli 

In this section, we focus on the evolution of the elastic properties of the material and neglect the 
plastic contribution to the granular behavior, that will be superimposed to the present analysis 
later in section [5l We hrst describe the numerical procedure to measure the elastic moduli of the 
anisotropic aggregate, and later we analyze the data and their relation with stress and fabric. 

4.1 Numerical probes 

In a general framework, a possible description for the incremental, elastic behavior of an 
anisotropic material is 


SP* 

where the isotropic and deviatoric components of stress have been isolated and are expressed as 
functions of and Sdev via a non-dimensional stiffness matrix |26| (by multiplying the moduli 



B 

Ai 


3fev 


A 2 

QOCt 


^'^dev 


( 11 ) 
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Figure 5: Evolution of (a) non-dimensional shear stress cr^gy and (b) deviatoric fabric Fdevj along 
the main strain path e^ev for the pure shear deformation mode, for volume fraction v = 0.706. 
The red symbols in (a - b) are the chosen states, which are first relaxed (blue ‘B’ symbols 
in (a) and (b)) and then used as initial configurations for the purely isotropic 2>5e^ and purely 
deviatoric fedev perturbations. 


with k*, the real stiffnesses can be extracted). B is the classical bulk modulus, and the 
octahedral shear modulus. The anisotropy moduli Ai and A 2 provide a cross coupling between 
the two parts (isotropic and deviatoric) of stress and strain increments. Ea. (|lll) provides a 
partial description for the evolving stress and stiffness of a sheared material, as it applies to 
a triaxial-box configuration (with eigensystem coincident with the axes of the box), where no 
shear strain/stress are measured and stress and moduli are assumed to be collinear. Moreover, 
the increase of stress and stiffness in the out-of-plane direction (z-direction here) due to the non- 
planar (triaxial) stress state associated with a the plane deformation mode, is not independently 
accounted for. These are rather hidden in the expression for deviatoric stress as proposed in 
Eq.Q and used in Ea. (|ll|) . However, we have chosen this representation, since advantages 
are obtained by investigating the elasticity of a granular material (e.g. soil), not through its 
resistance to direct stresses expressed by Young’s modulus and Poisson’s ratio, but rather in 
terms of (purely volumetric and deviatoric) stress-response to volume and shape changes, as 
described by the bulk modulus B and the octahedral shell modulus This aspect will be 

further addressed in section [5l where Ea. IjllD will be included in the theoretical model. 

To study the evolution of the effective moduli during shear, we choose different initial states 
(forty) as shown in Fig. [5l and apply sufficient relaxation, so that the granular assemblies 
dissipate the kinetic energy accumulated during the original shearing path. When the states 
along the shear path are relaxed, a much higher drop is visible in cr^g^ rather than in Fdev; see 
Fig. [5j This shows that the contact network remains almost intact and Fdev does not change; 
on the other hand, the average particle overlap is more sensitive to the relaxation stage and 
decreases, leading to a finite drop in cr^g^. Then we perform a small strain perturbation to 
these relaxed anisotropic states, i.e., we probe the samples, and measure the incremental stress 
response [50lll0]. Finally, the elastic moduli are calculated as the ratio between the measured 
increment in stress and the applied strain. We can obtain all the different moduli in Eq. 
applying an incremental pure volumetric or pure deviatoric strain and measuring the incremental 
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volumetric or shear stress response: 


B = 

^2 = 


5P* 


SSev 


SSSy 


^^dev—0 


^^dev—0 


= 


6 F* 


^^dev 


<5£v=0 


^oct _ ^'^dev 


^^dev 


<5ev=0 


( 12 ) 


Also for this part of the numerical experiment, the system is allowed to relax after the incremental 
strain is applied, that is the stress is measured after relaxation [5211501 Since the numerical probe 
experiments are carried out with zero contact friction, we are measuring the resistance of the 
frictionless material [40] , where only normal forces are involved. The first big question concerns 
the amplitude of the applied perturbation to get the elastic response [Tsmnii]. 


4.2 How small is small? 

In this section, we discuss the amplitude of the perturbations applied to measure the elastic 
stress response of the granular material. Also, we will discuss the results for larger amplitudes 
and the threshold between elastic and plastic regimes. 


4.2.1 Effect of isotropic perturbations 


Figs. [6] (column 1 and 2) show the changes in non-dimensional pressure SB*, non-dimensional 
shear stress isotropic fabric SFy and deviatoric fabric SFdev for different amplitudes of the 

isotropic perturbation applied to two relaxed states that have been sheared until ejev = 

0.0065 (nearly isotropic conhguration: column 1) and Sdev = 0.31 (steady state conhguration: 
column 2) respectively. The data correspond the the shear test with = 0.706. The linear 
elastic response is also plotted (red solid curve) in the whole strain range, as derived from the 
incremental behavior at very small strain, to give an idea of the deviation form elasticity when 
strain increases. 

SB* initially increases linearly and smoothly with SSEv, in agreement with the prediction of 
linear elasticity. Also the difference between the two initial states (near isotropic and steady 
state as shown in Figs. 6(a) and |6(b)' respectively) is minimal, meaning that the bulk modulus 
B (slope of SB* with SSEv in the elastic regime) is almost constant. This is not surprising, 
as we expect B to be dependent on isotropic quantities that, which stay mostly unchanged 
during the shear deformation, as discussed in section [T31 behaves similar as SB* for small 

strain, but shows several sharp drops for large strain. These correspond to sudden changes in 
the coordination number SC* (see Fig. CKa-b)), due to rearrangements in the system during 


the probe. For the nearly isotropic state (Fig. 6(e)), the ratio of (5crjgy with SSe^ in the linear 
elasticity regime, i.e. A 2 , is small when compared with the steady state (Fig. 6(f)). This clearly 


tells that A 2 evolves during the shear deformation for a given volume fraction, and must be 
linked with deviatoric quantities. 

SFy increases with SSSy, with more fluctuations compared to SB*, for both states considered 
here, Sdev = 0.0065 (nearly isotropic state, Fig. 6(i)) and Sdev = 0.31 (steady state. Fig. 6(j)). 
Moreover, the prediction using Eq. ([8]) for Fy, matches the dataset very well. Fdev does not 
change (<5Fdev = 0) with increasing SSsy, until the first rearrangement in structure occurs (see 
Figs. [7Kc-d)). After this <5Fdev starts to decrease with increasing amplitude 3<5ev, faster in the 
steady state (Fig. 6(m)) than in the near isotropic state, see Fig. 6(n), We note here that, 
when a non-incremental volumetric strain (SSSy > 10“^) is applied, the system moves from a 
volume-conserving to a new non-volume-conserving deformation path. As this system is already 
anisotropic, this leads to a decrease (SF^ev < 0) in deviatoric fabric Fdev, opposite to the increase 
{Sadev > 0) in deviatoric stress, see Figs. 6(e) and |6(f)[ higher in the steady state (Fig. 6(n)) 
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than in the nearly isotropic state (Fig. 6(m)). The last observation suggests that the distance 


between the volume conserving and non-volume conserving configurations increases with Sdev 
Hence, during isotropic compression (increasing of a pre-sheared (anisotropic) state, 

both the pressure P* and shear stress increase, with pressure increasing much faster leading 
to a decrease in deviatoric stress ratio Sdev = ^*devl^* ■ deviatoric fabric Fdev also decreases 
with isotropic compression of a pre-sheared state, and the decrease is initially faster than the 
exponential decay of Fdev (see section 0 below) with volume fraction v, as seen in Fig. 6(n) This 
decrease in Fdev becomes slower for large strain, as also seen in Fig. 6(m) These observations 
are consistent with the findings of Imole et al. |3Uj . where the authors noticed a decreasing 
steady state deviatoric fabric and deviatoric stress ratio with the increasing volume fraction, or 


4.2.2 Effect of deviatoric perturbations fede 


Figs.llKcolumn 3 and 4) show the changes in the same quantities as before for different amplitudes 
of the deviatoric perturbation <5edev) applied to a relaxed state with volume fraction v = 0.706 
that has been sheared until edev = 0.0065 (nearly isotropic configuration: column 3) and Sdev = 
0.31 (steady state configuration: column 4). 

6 P* increase linearly with Jsdev (the slope in the elastic regime is Hi), with Hi much smaller 


for the nearly isotropic state (Fig. 6(c)) than for the steady state (Fig. 6(d)). This shows that 
Hi evolves during the shear deformations, like H2, for a given volume fraction, and must be 
linked with the deviatoric state of the system. Moreover, after large deformation, both states 
show drops in 6 P*, which can be linked to the particle rearrangements at large deformation 
(see Fig. [7Kc-d)). A non-linear, irregular behavior shows up for Jedev > 10“^, with 6 P* positive 
in case of loose sample (present sample) and negative for dense samples (data not shown), in 
agreement with the observations in Fig. 2(a) also increases linearly with Jsdev) with 


(slope of the line) slightly higher for the near isotropic state (Fig. 6(g)) than for the steady state 


(Fig. 6(h)). Again, similar to 5P*, (^cr^g^ shows drops after large deformations, which can be 


linked to the particle rearrangements at large deformation (see Fig. [7Kc-d)). In the steady state, 
the incremental stresses {5P* and (^cr^g^) increase linearly for very small strain, as the relaxed 
configuration, starting point for the probes, has lower stress than the main deviatoric path (see 
Fig. |5(a)' ) and the system tends to regain the ’’missed” stress, when the shear restarts. After the 
first elastic response, SP* and (Jcr^g^ fluctuate around zero for larger amplitudes (Figs. 6(d)|and 
M ), as no change in stress is expected with increasing deviatoric strain in the steady state. 

6 Fv stays mostly zero when small dsdev is applied for both near isotropic and steady state 
configurations (Figs. 6(k)| and 6(1)). With increasing strain amplitude, 6 Fv increases in the case 
of a loose sample close to the isotropic state (Fig. 6(k)[), and decreases for denser samples (data 


not shown), in agreement with the behavior in Fig. 2(b) In Fig. 6(0), dP^ev for the nearly 
isotropic state, stays zero for <5edev < 10“^, when no rearrangements happen and the behavior is 
elastic, while it reaches a positive finite value for larger amplitude (that coincides with the slope 
of the curve in Fig. |3(b)[ ). This finite value increases with increasing anisotropy (or deviatoric 
strain state) until it reaches zero in the steady state, where no variation of deviatoric fabric is 
expected with further applied deviatoric strain (see Fig. |6(p)[ ). When compared to the model 
predictions in Ref. [30], the simulation data for Fdev well match with the theoretical line, where 
Fdev increases due to shear for the near isotropic state, and does not change for the steady state 
simulation. 
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4.2.3 Discussion and comparison 

Since we are interested in measuring the pure elastic response of the material, we take care that 
no rearrangements happen in the system during the numerical probe, that is 3(5ev and fedev are 
applied only up to 10“^ (with very slow wall movement rate ~ 10“®,i.e., smaller than for the 
main large shear strain preparation experiment). Looking at Fig. [6l we note that much bigger 
drops appear in the deviatoric response when the isotropic perturbation is applied. Vice-versa, 
the fluctuations/drops are much larger in pressure rather than in shear stress, when we deal 
with deviatoric perturbations. It is worthwhile to mention here that we have tested our method 
by applying strain perturbations in opposite directions i.e., 3<5ev and —36£v, Se^ev and —dedev 
This does not lead to any difference in the elastic response, as long as we stay in the limit of 
elastic perturbations. 

We test the rearrangements argument in Fig. [8l by plotting the calculated bulk modulus 
B and octahedral shear modulus against the amplitude of the applied isotropic 2>5£^ and 
deviatoric <5edev strain, respectively, for states at Edev = 0.0065 and 0.31 (nearly isotropic and 
steady state configurations, respectively) of the main deviatoric experiment. Both B and 
stay practically constant for small amplitudes and we can assume the regime to be linear elastic 
m- At 35 ev — 10 the first change in the number of contacts happens (Fig. EKa-b)) and B 
starts to increases non-linearly. Similarly, when Edev — 10“^, the first change in the number of 
contacts happens (Fig. EKc-d)) and starts to decay. It is interesting to notice that for both 
B and G°^^, the elastic regime shrinks when the main deviatoric strain Edev increases (Fig. [8]) 
and, also, when the volume fraction reduces, going towards the jamming volume fraction (data 
not shown). A similar modulus may be plotted for fabric as /Ssdev that, due to the finite 
size of the system, would be identically zero, until the first rearrangement occurs (see Fig. [^. 

We further check the elasticity of the probe by reversing the incremental strain. We plot the 
stress responses to volumetric/deviatoric strain in Fig. [9] and compare loading and unloading 
probes for different volume fractions (u = 0.706 and 0.812) and amplitudes. Looking at Figs. El 
[7] and [9] together, three regimes seem to appear. The first one for very small strain (< 5.10“®), 
due to the finite size of the system, is characterized by no opening and closing of contacts, and 
shows perfect reversibility of the data, i.e., elasticity in Figs. [9Ka-d). The second regime in 
Figs. llKe-h) shows some weakly irreversible behavior, but only for the smallest volume fraction 
and a mixed perturbation mode, see the sample at = 0.706 in Fig. |9(f)[ we associate this 
behavior to minor contact changes, as visible in Figs. [6] and El but no large scale rearrangements 
occur. Finally, the third regime, for perturbations two orders of magnitude higher (> 10“^), a 
residual strain after reversal shows up for both volume fractions and all types of perturbations, 
see Figs. EKi-l), proving also that plasticity is much more pronounced in the deviatoric modes 
than in isotropic ones. We claim that small drops are related to local (weak, almost reversible) 
re-structuring, while in the last case, the whole system (or big portion of it) is involved in the 
collapse of the structure, with a more pronounced effect for samples close to the jamming volume 
fraction [4^ 155] . 

For granular materials, the strain can not be split in elastic and plastic contributions by 
“trivially” referring to the residual deformation like in classical solids: as soon as we are out 
of the elastic range, rearrangements happen during loading and (even though less probably) 
during unloading, and most likely no original particle position is recovered. Finally, we note 
that the results shown here are valid for finite-size systems; for much larger (real) samples of 
much smaller particles, we expect the first elastic regime to reduce to much smaller strains. The 
boundary between the second and third regime is an issue for further research |68| . 
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4.3 Evolution of the moduli 

Using the four packings at different z/j, we next determine which variables affect the incremental 
response of the aggregates at different deviatoric strains along the main path. In order to 
understand the role of the microstructure, i.e., the fabric tensor F, the volumetric and deviatoric 
components, and F^ev, are considered. We postulate that the incremental response of the 
granular material can be uniquely predicted, once its fabric state (along with the stress state) 
is known, irrespective of the path that the system experienced to reach that state. In this sense 
the fabric tensor can be referred to as a state variable. 


4.3.1 Bulk modulus B 


In Fig. 10(a), we plot the incremental non-dimensional pressure 5P* against the amplitude 
of the applied isotropic perturbation 36ey for one volume fraction, u = 0.706, and various 
initial anisotropic configurations. The slope of each line is the bulk modulus of that state. It 
practically remains unchanged for different states and suggests that B is constant for a given 
volume fraction. 

In Fig. |10(b)[ we plot the variation of the bulk modulus B, with the isotropic fabric F^ for 
packings with different volume fractions t'j. B increases systematically when the five different 
reference configurations are compared, and it is related to the value of F^ constant at a given 
i^i ismoiiiD]. As expected B is a purely volumetric quantity and varies with changes in the 
isotropic contact network. The inset in Fig. |10(b)] shows that the bulk modulus remains almost 
constant with applied shear during a single deviatoric experiment m, behaving qualitatively 
similar to pressure P* and isotropic fabric Fy, see Figs. 1(a) and 3(a) respectively. That is, 
the contact orientation anisotropy, Fdevj which changes during the main deviatoric deformation 


path (see Fig. 3(b)) does not affect it. In agreement with observations on the volumetric fabric 


in section E21 also B shows a slight increase/decrease in the first part of the deviatoric path, 
more pronounced for loose samples, as clearly seen in Fig. |2(b)| The trend of B slightly deviates 
from TV in the low strain regime, while the dependence is well captured in the steady state, after 
large strain. The relation between bulk modulus and fabric was given in Ref. 


as: 


B = 


6 P* 


36e^ 


PoEv 


fedBv=0 . 


1 27 p ( Ev) + ( £v) (1 7p ( w)) 


91nFv 

<9(-ev) 


(13) 


where po, 'Jp and the jamming volume fraction i^c are fit parameters presented in Table [Tin 
pa 1.22 is dependent on the particle size distribution as presented in Refs. [MIESIEQ] , see 
section [2j For a given volume fraction, the above relation only requires the knowledge of the 
isotropic fabric F^ = g^i'C = g^i'C* (1 — (pr), where the empirical relations for C*{v) and </>r(z^) 
are taken from Refs. [301139] . see section [2j The numerical data show good agreement with the 
theoretical prediction presented in [23] and reported in Fig. 10(b) The minimum F^ is obtained 


at the jamming volume fraction, with = 0.658, C* = 6, and (j)r = ‘t’c = 0.13, leading to 
^mm _ ^ 2. At the jamming transition, we can extrapolate a finite value of the bulk modulus 
^mm _ Q 22, while it suddenly drops to zero below Uc [USl IS2 ESI ESI O EH EH ESj- The 
discontinuity of B is related to the discontinuity in F^, that jumps form zero to a finite value in 
I'c due to equilibrium requirements. 
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Modulus 

Fit parameter 

Bulk modulus B 

Po = 0.0425, 7p S3 0.2, Uc = 0.658 

First anisotropy modulus Ai 

ai = 1 ± 0.01 

Second anisotropy modulus A2 

ail = 1 i 0.02 

Octahedral shear modulus 

gi = 130 ± 3 


Table 1: Summary of fit parameters extracted from the small perturbation results in Eqs. m, 

(HH), disi), and ([I7|). 


4.3.2 Anisotropy moduli Ai and A 2 


In Fig. 11(a), we plot the non-dimensional pressure increment 6 P* against the strain amplitude, 
when the material is subjected to small deviatoric perturbations fedevi to measure the first 
anisotropy modulus Ai as defined in Eq. m, in given anisotropic configurations, as in Fig. 


10(a) Since the material is in an anisotropic state, an increment in deviatoric strain leads to 


a change in volumetric stress, along with shear stress. The slope of the curves, Ai, increases 
with the previous shear strain the system has experienced, going from small values in the initial 
isotropic configuration, to an asymptotic limit. 

We are interested in the ratio AijB. In this ratio, the dependence of isotropic fabric TV 
cancels out, all that remains is a pure dependence on Fjev In Fig. |ll(b)[ we plot the variation 
of AijB, with Tdev for packings with different volume fractions Ui as shown in the inset. Besides 
the fluctuations, the data collapse on a unique curve irrespective of volume fraction and pressure, 
that is, once a state has been achieved, a measurement of the overall anisotropy modulus is 
associated with a unique iViev An increasing trend of AijB with the fabric shows up. As the 


deviatoric fabric decreases with volume fraction (see Fig. 3(b)), this leads to lower values of the 


scaled anisotropy modulus for denser systems. In conclusion, we have a linear relation between 
for the first anisotropy modulus Aj: 

5P* 


Ai = 


— Oj\BF(\g 


<l£dev &^=0 
is the deviatoric part of fabric, and ai 


(14) 

0.66 is a fit parameter 


where B is the bulk modulus, Fde 
presented in Table [H 

we plot the stress response of the material to isotropic perturbation 

to measure the second 


12(a: 


In Fig. 

36£v, for the same anisotropic initial configurations as in Fig. 10(a) 
anisotropy modulus A2 as defined in Eq. m- Similarly to Ai, the slope of the elastic curves, 
i.e., A2, increases with the previous shear strain the system has felt, starting form zero until an 
asymptotic limit is reached. In Fig. 12(b), we plot the variation of A2/B, with F^ev for different 
volume fractions t'j as shown in the inset. Data show a very similar trend to what observed in 
Fig. |ll(b)| and besides the fluctuations, a collapse of data is observed]! Again we can relate A2 
to deviatoric fabric as: 

= aiiBFdev. (15) 

^^dev—0 


Ao = 


36 £v 


®Note that Uc for the same particulate system was reported as 0.66 for isotropic deformation in Ref. | 24| . 
as 0.6646 for isotropic and 0.658 for shear deformation in Ref. m- We use a similar i/c — 0.658 here, which, 
however, is dependent on history of the sample and on the deformation mode. The small deviations of B from 
Eq. (II. 8 II can be attributed to a (small) variation of Uc, however, this is beyond the focus of this paper. 

®A large data scatter is present in both figures Figs. 11(b) and 12(b)[ which increases for increasing deviatoric 
fabric fdev This is possibly due to other factors that may contribute to the evolution of the anisotropy moduli 
that are not considered in the present work. 
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The equality between the two fitting constants ai an ^ 1 (see Tabled]), states the symmetry 
of the stiffness matrix in Eq. (dU). 

Eq. dH and dSI provides an interesting, novel way to back-calculate the deviatoric structure 
in a granular sample via E^ev = where A and B can be inferred from wave propagation 

experiments, while the direct measurement of fabric is still an open issue [291EIIE11I36]. 


4.3.3 Octahedral shear modulus G' 


OCt 


In Eig. 13(a) 


tions in Eig. 


we plot the shear stress response of the material when the initial configura- 
10(a) are subjected to small deviatoric perturbations fedev The octahedral shear 
modulus is then measured, as defined in Eq. m- The slope of the curves for different 
initial configurations slightly decreases with the deviatoric state of the system, and saturates 
for high deformation edev; when the steady state is reached. Eig. |13(b)| shows the variation of 
^oct g^ga^i]2st shear strain edev starts from a finite value in the initial configuration, related 
to the isotropic contact network, and slightly decreases with increasing strain, with different 
rates for different volume fractions. The behavior of G°'^^ differs from that observed for the bulk 
modulus in the inset of Eig. |10(b) the shear resistance consistently decreases with shear strain 
and no transition between initial decrease/increase is observed, meaning that a factor other than 
Ey influences the change of G°^^ during the deviatoric path. Similarly to what done for Ai and 
A 2 , we look at the ratio of the shear modulus with the bulk modulus G°^^/B plotted against 


the isotropic fabric in Eig. 13(c), The ratio increases with increasing Ey, with higher values 


in the initial state than in the steady state (data are averaged over shear strain Sdev < 0.0065 to 
get the initial value and in the steady state to get the final one). The isotropic ratio (G'°‘^^/E). . 
increases with Ey, following the power law: 




1 — exp 


F _ pmni 

pet 


(16) 


where ^ ~ 0.51 represents the maximum value of ratio G°^^/B for large Ey (or volume 

fraction), E™'° ~ 4.2 is the volumetric fabric at the jamming transition, presented in section 
14.3.11 E" ~ 1.9 is the rate of growth of (G°^*/E).^., when the numerical data is extrapolated to 
the jamming transition, where (G°^*/E).^. = 0. This is in agreement with previous studies that 
find an upper limit equal 0.5 for the ratio between the shear and bulk moduli [3811731 fT8l I50| . 
In the limit of high Ey, the granular assembly becomes highly coordinated and practically 
follows the affine approximation that predicts a constant value for the ratio G°^^/B [78|. Here, 
a qualitatively similar behavior is observed for the values in the steady state, approaching a 
saturation ratio lower than the isotropic one. 

Next, in Eig. 13(d), we subtract the initial value (G°‘^*/E). . from G°^^/B and assume that 
Ey does not change during the deviatoric deformation. Interesting, we find that in this case the 
deviatoric microstructure alone is not able to capture the variation of the modulus along the 
shear path, but both stress a and fabric F seem to influence the incremental shear response, in 
agreement with findings in [H^. We relate the decrease of G°^^ to the deviatoric components of 
stress and fabric via: 


^oct _ ^^dev 


^^dev 


= B 


5£v=0 


Goct 


dl^deyBdev 


(17) 


where is the non-dimensional shear stress, Fdev is the deviatoric fabric and gi ~ 86 is a 
fit parameter reported in Table [T] Two contributions of the fabric to the shear stiffness can 
be recognized - isotropic and deviatoric. The overall contribution is multiplicative proportional 
to B, due to the isotropic contact network, changing very little with deviatoric strain. In the 
bracket, the first term gives the resistance of the material in the initial isotropic configuration. 
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whereas the second part only depends on the deviatoric (state) variables and characterizes 
the evolution of the shear modulus with deviatoric strain. That is, given the initial isotropic 
configuration, the corresponding is known |1611781 [50]; on the other hand, the deviation 
from isotropic to anisotropic network of such configuration uniquely defines the reduction in 
the shear stiffness. The joint invariant of deviatoric stress and fabric as proposed in 

izaiEzi, able to capture the evolution of the ratio of the elastic moduli along the whole undrained 
path, not only in the steady state, as seen in Fig. |13(d)[ 0 No more relation with volumetric 
quantities needs to be considered, as the evolution of cT^g.^,Tdev depends on the volume fraction 
of the sample Vi. 

a deviation from the fitting 


Note that when 0°^^^ is plotted against Eq. (flTli in Fig. 13(c) 


law is observed for each volume fraction, showing that extra correction terms might be needed 
for a more accurate description. This is neglected in this preliminary work. It is interesting to 
point out that the isotropic fabric has different effects in case of the anisotropy moduli Ai , A 2 
and as in the former two cases through B, is multiplied to Fdev and contributes to 

the growth of the moduli from zero to the asymptotic values, while in the latter case TV defines 
mostly the initial values of G°^^ via the bulk modulus, but does not affect the further decrease. 

In the next section, we use the evolution equation for the fabric as predicted from Eq. ([8]), 
and the relations between the elastic moduli and the stress and fabric, to predict an independent 
deformation experiment, namely the cyclic shear deformation, i.e., reverse shear after a large 
deviatoric strain. 


5 Prediction of undrained cyclic shear test 

In this section, the constitutive model is presented, involving the elastic moduli measured and 
calibrated in section S] and the plastic response of the material under large strain. The model 
is then used to predict the material response under cyclic shear, involving reversal. 


5.1 Calibration: Constitutive Model with Anisotropy 

We introduce here a constitutive model as proposed in Refs. [IH] EU ESI IMl EZj, extended to 
three dimensions, that takes into account the evolution of fabric, independently of stress: 


SP* = B36e^ + ASJede^, 

'^'^dev ~ A35£^ + G°‘^*5'CT(5e'dev) 

(5Fdev = /3Fsign(fedev)Ej()“S'i;'(5edev (18) 


In its simple, reduced form, the model involves only three moduli B, A and G°^^, defined in 
the previous section in Eqs. (nsj) - (nzi)- Due to A, the model provides a cross coupling between 
the two types of stress and strain in the model, namely the isotropic stress P* and shear stress 
cT^ev reacting to both isotropic (sy) and deviatoric (edev) strains. Fdev evolves differently from 
stress, as the rate of change with deviatoric strain can be (and in many cases is) different than 
the respective rate for the shear stress evolution. Note that additional terms (cross coupling of 
fabric with strain) might be needed for the incremental evolution of (5Edev in Eqs. (|18p . due to 
the observations from Fig. 6(n), where Edev and Ey change also with Sy and ffdev; respectively. 
However, both cross terms appear to be more activated in the highly anisotropic state, with 
values of the out-of-plane fabric considerably smaller than out-of-plane stress - but this has to 
be confirmed by other deformation paths also, i.e., we claim that some features are related to the 


^Such a split between isotropic and deviatoric fabric influence applies to this specific deformation path, where 
the volume is conserved. Additional terms may enter when non volume-conserving deformation paths are consid¬ 
ered. A very similar behavior is observed when the definition in Eq. ® is employed for the deviatoric fabric. 
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specific deformation path proposed here. If a dependence between stiffness and fabric similar to 
what proposed in Ea. ffT^ and (fT^ is assumed, previous arguments also lead to the conclusion 
that the out-of-plane stiffness terms developing during plane strain and neglected in Eo. dlip must 
be small compared io B, A and G. As a conclusion, for the sake of simplicity, both evolution 
of cross-coupling fabric terms and out-of-plane stiffness are neglected in the present work, and 
postponed to future investigations, for the description of arbitrary deformation paths. The use 
of non-frictional particles is another possible reason for the simplest model to work astonishingly 
well - so the general model is expected to show all contributions for arbitrary deformation, in 
the presence of friction. 

5,, = S/Si, with 5 = (1-Sdev/sS'eT) is a measure of the stress isotropy with normalized shear 
stress ratio Sdev = '^dev/-^*> is the initial stress isotropy at the start of a new deformation 

direction and/or after relaxation. 1 — 5o- is the measure for the probability of plastic events. 
Similarly, Sp = (1 — Tdev/-Tj^“)/5p is the fabric isotropy, and Sp is the initial fabric isotropy at 
the start of a new deformation direction and/or after relaxation. and represent the 

maximum (saturation) values of normalized shear stress ratio Sdev and deviatoric fabric Edev, 
respectively, and fip is the rate of change in Fdev at smaller strains (as shown in Eig. |3(b)[ ). 

It is worthwhile to point out that the definitions of S^ and Sp are different to those used in 
Refs. [48ll51j . as both 5^ and Sp are now scaled by the initial reference value and can take values 
between 0 and I. Due to S^j and Sp, the incremental response of the material is purely elastic, 
after relaxation or at strain reversal, with the elastic moduli evolving, as given by Eqs. m 
(ED, as functions of the momentary stress and structure states. At reversal, the probability for 
plastic deformation drops to zero and plastic events - as related to the approach to steady state 
- only occur after relatively large strain, that is the reversal stiffness is not affected. Due to 
Sa and Sp, the incremental response of the material in the large-strain steady state (S' = 0) 
becomes elastic (S = 1), just when the strain is reverted, or after relaxation (which is allowed 
before the probes). Due to the dependence of the elastic moduli on the stress/fabric state, 
the model involves non-linear elasticity in its present form (without contact non-linearities), 
while plasticity due to rearrangements is entirely associated to S^- On the other hand, the 
equation that describes the evolution of fabric is “purely plastic”, as there is no change in 
fabric (dE/ = 0), in the elastic regime, when no contact opening/closing and no multi-particle 
rearrangement happensU Thus the rate /3p is associated to changes of structure with deviatoric 
(shear) strain amplitude (not rate); changes are becoming more and more probably in the steady 
state. 

Now, we can predict an independent experiment, by using Eqs. (ED, and the relations for 
the four moduli B, A and with microscopic quantities given by Eqs. (ED “ (ED with the 
numerical scaling factors from Table [1] (starting from B, we can calculate the other moduli using 
the ratio). Moreover, four other parameters E)j^^ and /3p are needed to fully solve the 

coupled Eqs. (ED - The dependence of E)j^“ and I3p on volume fraction p, is well described 
by the exponential decay relation proposed in Refs. [301 EH], where constant values, as given in 
Fig. [TT]are used, as the volume is conserved during the cyclic shear test, as discussed next. 

®We want to point out here the difference between the non-linear elasticity built up along the main deviatoric 
path and the incremental elasticity, related to the small perturbations. Lets select two states A-B along the 
deviatoric path as indicated by points in Fig. [S] the incremental measured elastic response (moduli) is different 
between states A and B as it depends on stress and fabric, that is the stiffness matrix in Eq. ED , varies non- 
linearly with Edev On the other hand, when the incremental strain Ssdev is applied to each state (e.g., A or B), 
the incremental response is linearly elastic (by definition of incremental) and becomes plastic for high dsdev, as 
rearrangements happen and the moduli in that given state go from elastic to plastic. 
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5.2 Prediction: (Undrained) cyclic shear test 

We choose an initial isotropic configuration, with volume fraction v = 0.711 and apply deviatoric 
(volume conserving) shear for one cycle: loading, unloading and final re-loading, to recover 
the initial box configuration. Fig. [T3] shows the evolution of pressure P*, shear stress 
shear stress ratio Sdev and deviatoric fabric Pdev with deviatoric shear strain edev for one cycle, 
compared with the prediction using Eqs. m- Since the initial configuration is isotropic, the 
shear stress and Pdev start from zero and approach saturation values (with fluctuations) at 
large strains. During reversal, both drop with a soft response from their respective saturation 
value and decrease with unloading strain, crossing their zero values at different strain levels, 
and finally reach their steady state with negative signs. This supports the need of independent 
descriptions for the evolution of stress and fabric. Finally, re-loading is applied to reach the 
initial box configuration. 

The qualitative behavior of pressure P* is similar in simulations and model, going from a 
finite initial value to saturation with much less pronounced variations, since the deformation 
path is volume conserving. It is also interesting that the final state after the complete cycle, 
which corresponds to the initial box configuration, is highly anisotropic (non-zero stress 
and deviatoric fabric Pdev)- 

Both, the shear stress fj^g^ and deviatoric fabric Pdev, as well as their soft responses during 
strain reversal are well predicted by the model. P* increases during loading e^ev by ~ 9% and 
saturates at large strains. After reversal, P* drops because of opening and release of contacts and 
then increases again with unloading strain. Although P* is not quantitatively predicted by Eqs. 
m, the qualitative behavior is captured by the model, which requires a correction as proposed 
by Krijgsman and Ending m- The concept of a history dependent jamming point, introduced 
by Kumar et al. |TT], is capable of capturing the behavior of P* quantitatively, however, this 
goes beyond the scope of this study. 

Eqs. (fTHIl provide a set of equations able to describe the volumetric/deviatoric behavior of 
a granular assembly, in terms of stress and fabric. Once the initial state and the deformation 
path are defined, the evolution of isotropic fabric can be determined (using the coordination 
number and the fraction of rattlers) along the deformation path. The knowledge of isotropic 
and deviatoric fabric and the incremental relations in Eqs. (jlSp allow for the definitions of the 
moduli at each incremental step. Given also the probabilities for the plastic events (1 — (So- and 
1 — Sp), the coupled system can be solved. That is, the characterization of the initial state is the 
information needed to fully describe the behavior of the material along a general deformation 
path, defined in terms of strain, since the incremental evolution equations for both stress and 
structure are given. 

In the case of granular matter, the concept of a (homogeneous) material point in a continuum 
model is debated and many studies have been devoted to the introduction of a length scale in 
the constitutive model, starting from the Cosserat brothers, see miiMiiii] among others. Here 
we limit ourselves and state that a finite-size system is always needed, in order to calibrate any 
continuum model. That is, any model interpretation works only between the upper/lower bounds 
of infinite system and particle scale. When a finite-size system is considered an elastic range 
can always be detected, such that rearrangements happen (see section with negligible (tiny) 
probability for very small strain, and an elasto-plastic framework could then make sense. Here, 
we introduce a local rate-type model in Eqs. (|18p . and identify elasticity as the unique initial, 
static, configuration, from which the (incrementally irreversible) evolution of stress and structure 
follows. Our choice is to reduce elasticity to a “punctual range”, as plastic deformations (which 
include irreversible opening/closing of contacts by large scale rearrangements) will dominate for 
large deformations. Dynamics and kinetic fluctuations, leading to relaxation, are not considered 
here, but also needs to be taken into account, see e.g., [33] . 
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6 Summary and Outlook 


In a triaxial box, the four elastic moduli that describe the incremental, elastic constitutive 
behavior of an anisotropic granular material in terms of volumetric/deviatoric components, 
namely the bulk modulus B, the two anisotropic moduli Ai,A 2 and the octahedral shear modulus 
can be measured by applying small strain perturbations to relaxed states that previously 
experienced a large strain, volume conserving (undrained) shear path. A connection between the 
macroscopic elastic response and the micromechanics is established, by considering both stress 
and fabric tensors, a and F, respectively. While the bulk modulus B depends on the isotropic 
contact network FV, the deviatoric component of the fabric tensor Fdev is the fundamental state 
variable needed to properly model the ratios between the (cross-coupling) anisotropic and bulk 
moduli. When the deviatoric stress and strain are appropriately scaled (normalized), we find that 
the moduli reduce to three relevant ones, i.e. A = Ai = A 2 = aF^g^B. The anisotropy moduli 
are related to both deviatoric and isotropic fabric, as the whole contact network determines 
how the system will react to a perturbation. Surprisingly, when the shear resistance 
is considered, both the contact network and the deviatoric stress determine the incremental 
behavior of the assembly. When the initial response is subtracted, the residual ratio G°^^/B — 
(G°‘^*/i?).^. scales with the deviatoric state of the system, through the product cr^g^Tdev For 
strain amplitude larger than 10“^, rearrangements in the sample take place and the behavior 
deviates from elastic (reversible). The effect of increasing amplitude of isotropic/deviatoric 
strain perturbations on isotropic/deviatoric stress and fabric is investigated, in the case of nearly 
isotropic states and steady states at various different densities. For very small strain, the initial 
(linear) elastic regime, visible in the stress response, is associated to zero change in fabric. For 
higher strain amplitude applied to nearly isotropic state, plasticity comes into the play, and the 
incremental stress-strain relation deviates from linear as soon as the contact network changes. 
In the case of steady state, deviatoric strain can only induce fluctuations around the saturation 
value for both stress and fabric. Large volumetric strain induce substantial modifications, as 
the sample previously subjected only to volume-conserving deformation, experiences now large 
volume changes. In the limit of large strain, the tangential moduli of the stress-strain and 
fabric-strain curves (see Fig. [S]) are recovered. The relation between particle rearrangements 
and macro-scale plasticity is a present object of investigation, as well as the transition between 
local/global plastic regimes. As first important result, thanks to the independent study on 
elasticity, our study provides a new way to indirectly characterize the granular structure. Once 
the moduli in a given isotropic/anisotropic configuration, have been measured through wave 
propagation experiments, they can be uniquely associated with the internal fabric. However, we 
do not expect the proportionality to remain constant for different materials. 

As further step, a simple constitutive model is introduced that involves anisotropy, as pro¬ 
posed in Refs. |48ll51j . The non-linear elastic behavior is established and the irreversible/ plastic 
contribution is introduced via empirical probabilities for plastic events, that require more re¬ 
search and theoretical support. The dependence of the model parameters on volume fraction 
and polydispersity has been analyzed in previous extensive work [3D1 [3U]. Here, by using the 
new relations for the elastic moduli, we are able to integrate the increments at each state along 
a generic deformation path. Hence we can predict the evolution of pressure, shear stress and 
fabric for large strain, and also at and after reversal. The method is first calibrated and then 
applied to a volume conserving (undrained) shear cycle. When the prediction is compared with 
numerical simulations, quantitative agreement is found for the deviatoric field variables. The 
most notable feature of soft but different reversal responses of shear stress and fabric is well 
captured; the pressure response amplitude is underestimated by the present model. 

This study concerns a seemingly unrealistic material of spheres without friction and inter¬ 
acting with linear contact forces to exclude effects that are due to contact non-linearity, friction 
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and/or non-sphericity. This allows to unravel the peculiar interplay of stress with microstruc¬ 
ture. However, the work should be extended to more realistic cases involving particle shape, 
friction, and non-linear contact behavior. We expect that friction will not completely change 
qualitatively the observed relations between stiffness and fabric state, but possibly will add 
new effects to be explored in the future; the deviatoric fabric and the moduli are expected to 
change quantitatively when tangential forces are included. On the other hand, non-linearity at 
contacts will introduce an extra pressure-dependence for the moduli, as already shown by many 
authors (see e.g. [T6l ESllZl [50] in the case of Hertzian interactions). Speculating about the 
effects of shape goes beyond the scope of this study. A similar analysis is already in progress 
to check the influence of polydispersity on the relation between elastic stiffness and microstruc¬ 
ture, as polydispersity affects the contact network, the structure, and the orientation of contacts 

[22112511391- 

Future work will focus on the extension of our small perturbation approach to elasto- 
plasticity, by using concepts like e.g. the Gudehus response envelope [271152]. Other theoretical 
approaches involve ideas proposed by Einav or by Jiang and Liu [33], for which our results 
can provide a microscopically based calibration of parameters, but details are not discussed 
here. The information obtained for the pure elastic range can then be used to decouple the 
plastic contribution, associated with rearrangements, and to study the flow rule. The valida¬ 
tion of the present analysis with experimental data is another important goal. Nevertheless 
the issue of measuring fabric from laboratory experiments is far from solved, even though big 
advances have been made in recent years using photoelasticity, and microtomography CT-scans 
[291 [7211321 [5]. A partial validation is anyway possible when measuring the residual dependence 
of the elastic response from variables other than stress and porosity m, by means of acoustic 
measurements [36]. The behavior after more than one cycle deserves further investigation, from 
both simulational and theoretical points of view, to detect features like creep, liquefaction and 
ratcheting, analyzed in preliminary works m with constant elastic moduli and for many cycles 
m- Finally, a general tensor formulation that allows for highly different orientations of strain 
rate, stress and fabric is an open issue but can be inspired by the works of Thornton m and 
Zhao h Guo [57] . 
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Figure 6: (Rows) Change of non-dimensional pressure 6P*, non-dimensional shear stress ^cr^gy, 
isotropic fabric J-Fy and deviatoric fabric SF^ev versus strain amplitude. Column 1 and 2 repre¬ 
sent purely isotropic while column 3 and 4 represent deviatoric perturbation experiments. The 
perturbation is applied to the state corresponding to Edev = 0.0065 (nearly isotropic configura¬ 
tion; column 1 and 3) and Edev = 0.31 (steady state configuration: column 2 and 4) of the main 
deviatoric experiment with volume fraction u = 0.706. Note that the x-axis is log-scale, with 
inset plots in linear scale. The red line passing through the dataset in (a-j) represents a linear fit 
in the elastic regime for SSSy] de^ev < 10“^. The analytical predictions for the elastic range from 
our results section in Eqs. (fT3]l - (fT7|) are plotted as green line in (a-h). The green line in (i) 
and (j) represents Fy = gsi^C calculated using Eq. ([8]), when subtracted from its initial value. 
The dashed horizontal line in (k)-(p) represents zero. The green line in (m) and (n) represent 
the evolution of change in deviatoric fabric <5Fdev in critical state using parameters from Table 
3 of Ref. [30], with the assumption that the new state after volumetric deformation is also in 
critical state. The green line in (o) and (p) represents Eq. (18) from Ref. [3U] when subtracted 
from its initial value F^g^ = 0.03 for (o) and = 0.113 for (p), with the growth rate Pf = 39 
and = 0.12. 
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Figure 7: Change of the coordination number 5C* = <5 {M^/N^) (black curve) and the 
modified coordination number 6C* = 5{M^/Ni) (red curve), defined in section [2l versus 
strain amplitude during purely (a-b) isotropic and (c-d) deviatoric perturbation experiments 
(corresponding plots as in Fig. [U]). The perturbation is applied to the state corresponding 
to edev = 0.0065 (nearly isotropic configuration: (a) and (c)) and Sdev = 0.31 (steady state 
configuration: (b) and (d)) of the main deviatoric experiment with volume fraction v = 0.706. 
Note that the x-axis is on log-scale, with inset plots in linear scale. 



Figure 8: Evolution of (a) bnlk modulus B and (b) octahedral shear modnlus with the 
respective applied isotropic Sfev and deviatoric fedev strain amplitudes for a state corresponding 
to edev = 0.0065 (nearly isotropic configuration: green ‘B’) and edev = 0.31 (steady state 
configuration: blue ‘•’) of the main deviatoric experiment with volume fraction u = 0.706. 
Corresponding dashed horizontal lines represents the initial values of B and G°^^. 
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Figure 9: Evolution of non-dimensional pressure P*, non-dimensional shear stress during 
small (a - d), medium (e-h), and large (i - 1) perturbations in the loading (symbols) and 
then unloading (solid lines) direction. Red ‘-I-’ represents loading and the green line represents 
unloading for v = 0.812. Similarly, blue represents loading and the black line represents 
unloading for v = 0.706. The deformation is applied to the state corresponding to ejev = 0.31 
(steady state configuration) of the main deviatoric experiment. 
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Figure 10: (a) Evolution of change in non-dimensional pressure 6P* during purely isotropic 

perturbations Shev for different states for volume fraction u = 0.706 along the main path as 
shown in the inset, (b) Evolution of the bulk modulus B as scaled with isotropic fabric for 
five different volume fraction as shown in the inset. The solid line passing through the data 
represents Eq. m- The dashed lines represent the initial and steady state data, as given in the 
legend. 
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Figure 11: (a) Evolution of change in non-dimensional pressure 6P* during purely deviatoric 

perturbations (Jffdev for different states for volume fraction = 0.706 along the main path as 
shown in the inset. The arrow indicates the direction of increasing strain states during main 
deviatoric experiments, (b) Evolution of the ratio of first anisotropy modulus with bulk modulus 
AijB as function of the deviatoric fabric F^ev for five different volume fractions as shown in the 
inset. The solid line passing through the data represents Eq. (fTT)l divided by B. 
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Figure 12: (a) Evolution of change in non-dimensional shear stress during purely isotropic 
perturbations 3fev for different states for volume fraction v = 0.706 along the main path as 
shown in the inset. The arrow indicates the direction of increasing strains during main deviatoric 
experiments, (b) Evolution of the ratio of second anisotropy modulus with bulk modulus A 2 IB 
as scaled with the deviatoric fabric Fjev for five different volume fraction as shown in the inset. 
The solid line passing through the data represents Eq. (USD divided by B. 
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Figure 13: (a) Change in shear stress versus strain amplitude during purely deviatoric 

perturbations fedev for different states, with volume fraction v = 0.706, along the main path as 
shown in the inset, (b) Evolution of octahedral shear modulus along the main deviatoric 
path edev for five different volume fractions as shown in the inset. The corresponding lines 
passing through the data represents Eq. (HZl). (c) Evolution of ratio of octahedral shear modulus 
and bulk modulus, i.e., G°^^/B with isotropic fabric F^, together with the averaged values at 
the initial (near isotropic state averaged over shear strain Sdev < 0.0065) and the steady state 
(averaged dataset in the steady state), as given in the legend. Note that the difference between 
initial and steady state increases with denser systems. The solid orange line passing through 
the isotropic dataset represents Eq. (1161) . (d) Evolution of the ratio of octahedral shear modulus 
and bulk modulus when its initial value, i.e., G°^^/B — {G°^^/B^ . is subtracted, plotted using 
Eq. (I17p . for five different volume fractions as shown in the inset. 
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Figure 14: Evolution of (a) pressure P*, (b) stress (c) normalized stress Sdev, and (d) 

deviatoric fabric Fdev with shear strain ejev during cyclic shear at constant volume i' = 0.711, 
starting from an initial isotropic configuration. The values of and Pf for i' = 0.711 

are 0.167, 0.124 and 40.04 respectively, taken directly from the relations proposed in Refs. 
1301 [39]. The red data points are the DEM simulation data over which the calibration of 
moduli was done, while the green data points represents unloading (reversal) and re-loading. 
The solid line is the prediction of the DEM observations using Eqs. m- 
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